#include <stdlib.h>
#include <stdio.h>
#include "nbody.h"
#include "video.h"
#include <time.h>
#include <math.h>
#include <string.h>
#include <sys/time.h>
#include "mpi.h"

double *x, *y; 
double *fx, *fy, *vx, *vy, *m;
int    nbodies;
int    width;
int    height;
int    nsteps;
double g;
double interval;
double min_dist;

int    **grid;

/**
 *	Inicializa los datos criticos de los cuerpos existentes en el fichero de entrada
 */
int init_bodies(char *file) {
	FILE *f; 
	int i;
 
	f = fopen(file,"r");
	
	fscanf(f,"%d %lf %lf %d %lf %d %d\n",&nbodies,&g,&interval,&nsteps,&min_dist,&width,&height);
	x  = malloc(sizeof(double) * nbodies);
	y  = malloc(sizeof(double) * nbodies);
	fx = malloc(sizeof(double) * nbodies);
	fy = malloc(sizeof(double) * nbodies);
	vx = malloc(sizeof(double) * nbodies);
	vy = malloc(sizeof(double) * nbodies);
	m  = malloc(sizeof(double) * nbodies);
	for (i = 0; i < nbodies; i++) {
		fscanf(f,"%lf %lf %lf %lf %lf\n",&x[i],&y[i],&vx[i],&vy[i],&m[i]);
	}
	fclose(f);
	printf ("NBODY PROBLEM____________________________________________________________________\n");
	printf ("nbodies: %d\tg var: %0.2lf\tinterval: %0.2lfsec\tnsteps: %d\tmin_dist: %0.2lfp\twidth %d\theight: %d\n", nbodies,g,interval,nsteps,min_dist,width,height);
	return 1;
}

/**
 *	Escribe el fichero de salida tras realizar todos los calculos
 */
int finalize_nbodies(char *file) {
	FILE *f;
	int i;
     
	f = fopen(file,"w");
	fprintf(f,"%d %lf %lf %d %lf %d %d\n",nbodies,g,interval, nsteps,min_dist,width,height);
	for (i = 0; i < nbodies; i++) {
		fprintf(f,"%lf %lf %lf %lf %lf\n",x[i],y[i],vx[i],vy[i],m[i]);
	}
	fclose(f);
	free(x);
	free(y);
	free(fx);
	free(fy);
	free(vx);
	free(vy);
	free(m);
	return 1;
}

double factorial(double x) {
	double fac = 1;
	int i;
	for (i=2; i<=x; i++) {
		fac *= i;
	}
	return fac;
}

double mypow(double x, double exp) {
	double result = 1;
	int i;

	result = 1;
	for (i=0; i<exp; i++) {
		result *= x;
	}
	
	return result;
}

double mysin(double alpha) {
	double result = 0;
	int i;

	for (i=0; i < MAX_ACCURACY;i++) {
		result += (double)(mypow(-1,i) / factorial((double)(2*i+1))) * (double)mypow(alpha, 2*i+1);
	}
	return result;
}

double mycos(double alpha) {
	double result = 0;
	int i;

	for (i=0; i < MAX_ACCURACY;i++) {
		result += (double)(mypow(-1,i) / factorial((double)(2*i))) * (double)mypow(alpha, 2*i);
	}
	return result;
}


int calc_nbodies() {
	double *ax, *ay;
	double fax, fay;
	double dist, f, alpha;
	int i,j;     

	ax = malloc(sizeof(double) * nbodies);
	ay = malloc(sizeof(double) * nbodies);

	//Pone la diagonal a cero
	for (i = 0; i < nbodies; i++) {
		fx[i] = fy[i] = 0;
	}
	
	//Debemos tener en cuenta que nos haya llegado todos los mensajes de MPI
	for (i = 0; i < nbodies; i++) {
		for (j = i + 1; j < nbodies; j++) {
			dist = sqrt((x[j] - x[i]) * (x[j] - x[i]) + (y[j] - y[i]) * (y[j] - y[i]));
			if (dist > min_dist ) {
				
				//No tocar todo esto
				f = g * m[i] * m[j] / (dist * dist);

				if (f > MAX_F) {
					f = MAX_F;
				}
	
				if (y[j]-y[i] == 0) {
			
					alpha = atan(0);
			
				} else if (x[j]-x[i] == 0) {
			
					if (y[j]-y[i] > 0) {
						alpha = M_PI_2;
					} else {
						alpha = - M_PI_2;
					}
			
				} else {
			
					alpha = atan(y[j]-y[i] / x[j]-x[i]);
			
				}
		
				fax   = f * mycos(alpha);
				fay   = f * mysin(alpha);
				fx[i] += fax;
				fy[i] += fay;
				fx[j] -= fax; 
				fy[j] -= fay;
			}
		}
	}
	for (i = 0; i < nbodies; i++) {
		ax[i] = fx[i] / m[i];
		ay[i] = fy[i] / m[i];
		vx[i] = vx[i] + ax[i] * interval;
		vy[i] = vy[i] + ay[i] * interval;
		x[i] = x[i] + vx[i] * interval;
		y[i] = y[i] + vy[i] * interval;
		 
		if (x[i] <= 0) {
			x[i] = 2;
			vx[i] = - vx[i];
		}
		if (x[i] >= width ) {
			x[i] = width - 2;
			vx[i] = - vx[i];
		}
		if (y[i] <= 0) {
			y[i] = 2;
			vy[i] = - vy[i];
		}
		if (y[i] >= height) {
			y[i] = height - 2;
			vy[i] = - vy[i];
		}
		//printf("body %d: x=%d, y=%d vx=%lf vy=%lf ax=%lf ay=%lf fx=%lf  fy=%lf\n", i, x[i], y[i],vx[i],vy[i],ax[i],ay[i],fx[i],fy[i]);
	}

	free(ax);
	free(ay);

	return 1;
}



//int calc_nbodies_nuevo(int cuerpo, int columna) {
int calc_nbodies_nuevo(char *cuerpos) {
	//double *ax, *ay;
	//double fax, fay;
	double /*dist,*/ f, alpha;
	//int i,j;     

	double distancia = 0;
	
	//ax = malloc(sizeof(double) * nbodies);
	//ay = malloc(sizeof(double) * nbodies);

	int xi=0, xj=0, yi=0, yj=0, mx=0, my=0;
	
	printf("CUERPOS\n%s\n\n\n", cuerpos);
	
	
	printf("*******probando*****\n");
	
	sscanf(cuerpos, "%d %d %d %d %d %d %lf", &xi, &yi, &mx, &xj, &yj, &my);
	
	printf("_______ xi: %d\n", xi);
	printf("_______ xj: %d\n", xj);
	printf("_______ yi: %d\n", yi);
	printf("_______ yj: %d\n", yj);
	
	distancia = sqrt((xj - xi) * (xj - xi) + (yj - yi) * (yj - yi));
	
	printf("Distancia entre los cuerpos (%d,%d) y (%d,%d) es: %lf\n", xi, yi, xj, yj, distancia);
	printf("Distancia minima: %lf\n", min_dist);
	
	
	if (distancia > min_dist) {
		f = g * mx * my / (distancia * distancia);
		printf("Fuerzas entre cuerpos: %lf\n", f);
	}
	
	
	
	
	
	//Debemos tener en cuenta que nos haya llegado todos los mensajes de MPI
	/*for (i = 0; i < nbodies; i++) {
		for (j = i + 1; j < nbodies; j++) {
			dist = sqrt((x[j] - x[i]) * (x[j] - x[i]) + (y[j] - y[i]) * (y[j] - y[i]));
			if (dist > min_dist ) {
				
				//No tocar todo esto
				f = g * m[i] * m[j] / (dist * dist);

				if (f > MAX_F) {
					f = MAX_F;
				}
	
				if (y[j]-y[i] == 0) {
			
					alpha = atan(0);
			
				} else if (x[j]-x[i] == 0) {
			
					if (y[j]-y[i] > 0) {
						alpha = M_PI_2;
					} else {
						alpha = - M_PI_2;
					}
			
				} else {
			
					alpha = atan(y[j]-y[i] / x[j]-x[i]);
			
				}
		
				fax   = f * mycos(alpha);
				fay   = f * mysin(alpha);
				fx[i] += fax;
				fy[i] += fay;
				fx[j] -= fax; 
				fy[j] -= fay;
			}
		}
	}
	for (i = 0; i < nbodies; i++) {
		ax[i] = fx[i] / m[i];
		ay[i] = fy[i] / m[i];
		vx[i] = vx[i] + ax[i] * interval;
		vy[i] = vy[i] + ay[i] * interval;
		x[i] = x[i] + vx[i] * interval;
		y[i] = y[i] + vy[i] * interval;
		 
		if (x[i] <= 0) {
			x[i] = 2;
			vx[i] = - vx[i];
		}
		if (x[i] >= width ) {
			x[i] = width - 2;
			vx[i] = - vx[i];
		}
		if (y[i] <= 0) {
			y[i] = 2;
			vy[i] = - vy[i];
		}
		if (y[i] >= height) {
			y[i] = height - 2;
			vy[i] = - vy[i];
		}
		//printf("body %d: x=%d, y=%d vx=%lf vy=%lf ax=%lf ay=%lf fx=%lf  fy=%lf\n", i, x[i], y[i],vx[i],vy[i],ax[i],ay[i],fx[i],fy[i]);
	}

	free(ax);
	free(ay);
	*/

	return 1;
}




 
int main (int argc, char ** argv) {
	int iter;
	int x11 = 0;
	int tipo_reparto = 0;
	//struct timeval start, end;
	double usec = 0;
	
	/**
	 *
	 *
	 */
	int rango;
	int numero_procesos;
	int numero_interacciones = 0;
	int grano = 0;
	int resto = 0;
	int contadorCuerpos = 0;
	int contadorProcesos = 0;
	int contadorInteracciones = 0;
	int contadorGrano = 0;
	int contadorColumnas = 0;
	int hay_trabajo = 1;


	int *env;

	
	
	char *enviados;
	char *recibidos;
	
	
	
	MPI_Init(&argc, &argv);
	MPI_Comm_rank(MPI_COMM_WORLD, &rango);
	MPI_Comm_size(MPI_COMM_WORLD, &numero_procesos);
	MPI_Status status;
	

	//iter=0;
	if (rango == 0){
		if (argc < 5) {
			printf ("Use:\n");
			printf ("   %s <fichero_entrada> <fichero_salida> <X11_0|1> <tipo_reparto>\n", argv[0]);
			printf ("Donde:\n");
			printf ("   <X11_0|1> pueder ser o bien '0' (ejecuci�n en modo texto) o bien '1' (ejecuci�n en modo gr�fico).\n");
			printf ("   <tipo_reparto> pueder ser o bien '0' (reparto balanceado) o bien '1' (reparto no balanceado).\n");
			return 0;
		}

		// Comprobamos valor X11 dado por parametro
		x11 = atoi(argv[3]);
		// Si es menor que cero o mayor que uno, entonces fin del programa
		if(x11 < 0 || x11 > 1) {
			printf("ERROR: El valor X11 solo puede ser '0' (modo texto) o '1' (modo gr�fico).\n");
			return 1;
		}

		// Comprobamos valor tipo_reparto pasado por parametro
		tipo_reparto = atoi(argv[4]);
		// Si es menor que cero o mayor que uno, entonces fin del programa
		if (tipo_reparto < 0 || tipo_reparto > 1){
			printf("ERROR: El valor del tipo de reparto solo puede ser '0' (reparto balanceado) o '1' (reparto no balanceado).\n");
			return 1;
		}

		//inicializamos las variables criticas en funci�n del fichero de entrada
		init_bodies(argv[1]);

		//Si es particionamiento balanceado calculamos el numero de cuerpos para cada proceso
		if (tipo_reparto == 0){
			//Calculamos numero de interacciones desde Cuerpo 1 hasta Cuerpo N
			int j;

			for (j = 1; j < nbodies; j++){
				numero_interacciones = numero_interacciones + nbodies - j;
			}
			
			printf("Numero interaciones: %d\n", numero_interacciones);

			grano = numero_interacciones/numero_procesos;
			
			//Si la divisi�n es exacta...
			if ((numero_interacciones%numero_procesos) == 0){
				printf("Division exacta del grano....\n");
				resto = 0;
			} else { //si no es exacta
				printf("Mas uno al grano....\n");
				
				resto = numero_interacciones%numero_procesos;
			}
			

			printf("Grano = %d\nResto = %d\n", grano, resto);
			printf("Nbodies = %d\n", nbodies);
			printf("Num_proc = %d\n", numero_procesos);
			
		} else {
			
		}
		
		//Reservamos memoria con tamaño igual al número de cuerpos por proceso junto con los cuerpos que puedan sobrar
		
		env = malloc((numero_interacciones * (grano + resto)) * sizeof(int));
		
		
		
		//enviados = malloc((numero_interacciones * (grano + resto)) * sizeof(char));
		enviados = malloc((numero_interacciones * (nbodies * (grano + resto))) * sizeof(char));
		//recibidos = malloc((numero_interacciones * (grano + resto)) * sizeof(char));
		recibidos = malloc((numero_interacciones * (nbodies * (grano + resto))) * sizeof(char));
		
		
	}
	
	MPI_Bcast (&nbodies, 1, MPI_INT, 0, MPI_COMM_WORLD);
	MPI_Bcast (&grano, 1, MPI_INT, 0, MPI_COMM_WORLD);
	MPI_Bcast (&resto, 1, MPI_INT, 0, MPI_COMM_WORLD);
	MPI_Bcast (&numero_interacciones, 1, MPI_INT, 0, MPI_COMM_WORLD);
	MPI_Bcast (&min_dist, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
	MPI_Bcast (&g, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
	MPI_Bcast (&interval, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
	MPI_Bcast (&hay_trabajo, 1, MPI_INT, 0, MPI_COMM_WORLD);
	MPI_Barrier(MPI_COMM_WORLD);	
	

	//Se activa modo gr�fico
	if (x11 == 1) {
		allocate_grid(width, height);
		setupWindow(width, height);
	}
	
	
	
	if (rango == 0) {
		
		if (tipo_reparto == 0) {
			
			for (iter = 0; iter < nsteps; iter++) {
				//Para saber cuantos cuerpos totales hay
				contadorCuerpos = 0;
				//Para saber cuanto procesos hay
				contadorProcesos = 0;
				//Para saber cuantas iteraciones hay
				contadorInteracciones = 0;
				//Para saber cuantos cuerpos se han pasado ya
				contadorGrano = grano;
				
				contadorColumnas = 0;
				
				//Mientras no esten todos los cuerpos tratados....
				while (contadorCuerpos < (nbodies-1)) {
					//Para cada proceso....
					for(contadorProcesos = 0; contadorProcesos < numero_procesos; contadorProcesos++) {
						//Si proceso 0
						if (contadorProcesos == 0){
							//el resto va para el primer proceso
							contadorGrano += resto;
							//Contabilizamos hasta tener el grano completo
							while (contadorGrano > 0) {
								//Recorremos la fila hasta llegar a la última columna
								if (contadorColumnas > contadorCuerpos) {
									//almacenamos la posicion que le corresponde
									
									strcpy(recibidos, "");
									sprintf(recibidos, "%d %d %d %d %d %d", (int)x[contadorCuerpos], (int)y[contadorCuerpos], (int)m[contadorCuerpos], (int)x[contadorColumnas], (int)y[contadorColumnas], (int)m[contadorColumnas]);
									printf("RECIBIDOS\n%s\n\n\n", recibidos);
									
									
									printf("proceso 0 hace calc_nbodies_nuevo.....................\n");
									calc_nbodies_nuevo(recibidos);
									printf("proceso 0 termina calc_nbodies_nuevo.....................\n");
									
									//ya nos queda menos para completar el grano completo
									contadorGrano--;
									//ya hemos realizado una interaccion mas de las totales
									contadorInteracciones++;
									//si hemos hecho todas las interacciones
									if (contadorInteracciones > (numero_interacciones-1)) {
										//paramos
										hay_trabajo = 0;
									}
									
								}
								//Pasamos a una nueva columna
								contadorColumnas++;
								//Si hemos recorrido todas las columnas
								if (contadorColumnas >= nbodies) {
									//Nueva fila
									if (contadorCuerpos >= nbodies){
										hay_trabajo = 0;
									} else {
										contadorCuerpos++;
									}
									//Reiniciamos columnas
									contadorColumnas = 0;
								}
							}
							//Reiniciamos el grano
							contadorGrano = grano;
						//Resto de procesos
						} else {
							//Le corresponde un grano
							while (contadorGrano > 0) {
								//Recorremos la fila hasta llegar a la última columna
								if (contadorColumnas > contadorCuerpos) {
									//almacenamos la posicion que le corresponde									
									
									strcpy(enviados, "");
									sprintf(enviados, "%d %d %d %d %d %d", (int)x[contadorCuerpos], (int)y[contadorCuerpos], (int)m[contadorCuerpos], (int)x[contadorColumnas], (int)y[contadorColumnas], (int)m[contadorColumnas]);
									printf("ENVIADOS a P%d\n%s\n\n\n", contadorProcesos, enviados);
									
									MPI_Send(&enviados, numero_interacciones * (grano + resto), MPI_CHAR, contadorProcesos, 1, MPI_COMM_WORLD);
									printf("enviado correctamente......_____ %s  ______-..........\n", enviados);
									
									//ya nos queda menos para completar el grano completo
									contadorGrano--;
									//ya hemos realizado una interaccion mas de las totales
									contadorInteracciones++;
									//si hemos hecho todas las interacciones
									if (contadorInteracciones > (numero_interacciones-1)) {
										//paramos
										hay_trabajo = 0;
									}
									
								}
								//Pasamos a una nueva columna
								contadorColumnas++;
								//Si hemos recorrido todas las columnas
								if (contadorColumnas >= nbodies) {
									//Nueva fila
									if (contadorCuerpos >= nbodies){
										hay_trabajo = 0;
									} else {
										contadorCuerpos++;
									}
									//Reiniciamos columnas
									contadorColumnas = 0;
								}
							}
							//Reiniciamos el grano
							contadorGrano = grano;
						}
						
						//printf("%s\n\n\n", enviados);
						
						/*if (contadorProcesos > 0) {
							//MPI_Send(enviados, numero_interacciones * (grano + resto), MPI_CHAR, contadorProcesos, 1, MPI_COMM_WORLD);
							printf("Envio %d datos\n", numero_interacciones *(grano + resto));
							
							MPI_Send(&contadorCuerpos, numero_interacciones * (grano + resto), MPI_INT, contadorProcesos, 1, MPI_COMM_WORLD);
							MPI_Send(&contadorColumnas, numero_interacciones * (grano + resto), MPI_INT, contadorProcesos, 1, MPI_COMM_WORLD);
							
							
							//MPI_Send(&posicion[contadorCuerpos][contadorColumnas], numero_interacciones * (grano + resto), MPI_INT, contadorProcesos, 1, MPI_COMM_WORLD);
							//MPI_Send(&x, numero_interacciones * (grano + resto), MPI_INT, contadorProcesos, 1, MPI_COMM_WORLD);
							
							printf("Enviado el cuerpo %d y %d\n", contadorCuerpos, contadorColumnas);
							
						}*/
						
					}
					
					//printf("LLEGO\n");
					/*for (i = 1; i < numero_procesos; i++){
						//printf("ENTRO\n");
						//int *recv = malloc((numero_interacciones * (grano + resto)) * sizeof(int));
						int recv;
						MPI_Recv((int *)recv, (grano + resto) * numero_interacciones, MPI_INT, i, 1, MPI_COMM_WORLD, &status);
						printf("P%d->Mensaje_Recibido de %d:\t%d||\n",rango, i, recv);
					}
					printf("SALGO\n");*/
					
				}
				
				
				
			}
			/**for (iter = 0; iter < nsteps; iter++) {
				//gettimeofday(&start, NULL); 
				calc_nbodies();
				//gettimeofday(&end, NULL); 
				//usec += end.tv_usec + 1000000*(end.tv_sec - start.tv_sec) - start.tv_usec;
				if (x11) {

					init_grid(width, height);

					for (i=0; i < nbodies; i++) {
						grid[(int)x[i]][(int)y[i]]= m[i];
					}

					if(iter % REFRESH ==0 ) {
						myDraw(width, height);
					}
				}
			}*/
			
		} else {
			
		}
	} else {
		printf("Voy a recibir....\n");
		
		printf("----------- hay_trabajo = %d -------------\n", hay_trabajo);
		
		
		while (hay_trabajo != 0) {
			printf("P%d->Esperando\n", rango);
			printf("Recivo %d datos\n", numero_interacciones * (nbodies * (grano + resto)));
			
			char *rcv = malloc((numero_interacciones * (nbodies * (grano + resto))) * sizeof(char));
			
			printf("numero_interacciones ==== %d\n", numero_interacciones);
			printf("nbodies ==== %d\n", nbodies);
			printf("grano ==== %d\n", grano);
			printf("resto ==== %d\n", resto);
			
			printf("^^^^^^^^^^^^^\n");
			printf("_-¬¬¬¬¬¬ RECIBIDO ¬¬¬¬¬¬-_ %s\n", recibidos);
			printf("^^^^^^^^^^^^^\n");
			
			//MPI_Recv(&recibidos, numero_interacciones * (nbodies * (grano + resto)), MPI_CHAR, 0, 1, MPI_COMM_WORLD, &status);
			MPI_Recv(&rcv, numero_interacciones * (nbodies * (grano + resto)), MPI_CHAR, 0, 1, MPI_COMM_WORLD, &status);
			
			printf("RECIBIDO ==== %s\n", rcv);
			
			
			printf("proceso %d ===== calculando calc_nbodies_nuevo.......\n", rango);			
			
			calc_nbodies_nuevo(recibidos);
			
			printf("proceso %d ===== calc_nbodies_nuevo calculado.......\n", rango);

		}
	}

	if (x11 == 1) {
		cleanup_grid(width, height);
	}

	//Sincronizacion de procesos antes de liberar recursos
	MPI_Barrier(MPI_COMM_WORLD);
	//Proceso 0 libera recursos
	if (rango == 0){
		finalize_nbodies(argv[2]);
		free(enviados);
		free(recibidos);
		printf ("nbody takes: %0.2fsec    mean: %0.2fsec    accuracy: %d \n", (double)usec / 1000000, (double)usec / 1000000/nsteps,MAX_ACCURACY);
	}

	MPI_Finalize();
	
	return 1;
}
